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ABSTRACT 

The flow resulting from a strong ultrarelativistic shock moving through a stel- 
lar envelope with a polytrope-like density profile has been studied analytically 
and numerically at early times while the fluid temperature is relativistic — that is, 
just before and just after the shock breaks out of the star. Such a flow should ex- 
pand and accelerate as its internal energy is converted to bulk kinetic energy; at 
late enough times, the assumption of relativistic temperatures becomes invalid. 
Here we present a new self-similar solution for the post-breakout flow when the 
accelerating fluid has bulk kinetic Lorentz factors much larger than unity but is 
cooling through p/n of order unity to subrelativistic temperatures. This solution 
gives a relation between a fluid element's terminal Lorentz factor and that ele- 
ment's Lorentz factor just after it is shocked. Our numerical integrations agree 
well with the solution. While our solution assumes a planar flow, we show that 
corrections due to spherical geometry are important only for extremely fast ejecta 
originating in a region very close to the stellar surface. This region grows if the 
shock becomes relativistic deeper in the star. 



1. Introduction 

The energy and Lorentz factor that we expect in the ejecta in supernovae and gamma- 
ray bursts are important because they constrain the amount of energy that can be deposited 
in the photons we observe from these explosions. Previous work on the ejecta, notably Tan 
et al. (2001), uses as a starting point the analytic solutions of Johnson and McKee (1971) 
for a planar relativistic shock propagating into cold surroundings: by the time the shock 
reaches the outer envelope of the star, the likely source of the ejecta, it has accelerated to 
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relativistic speeds and its geometry is planar. Several authors have described the reltivistic 
shock's propagation and acceleration through this envelope (see, for example, Best and Sari 
2000; Perna and Vietri 2002; Sari 2006), but the work of Johnson and McKee (1971) and 
other analytic work on the flow from a relativistic shock that breaks out of a star (Nakayama 
and Shigeyama 2005; Pan and Sari 2006) show that significant acceleration also occurs after 
the fluid is shocked. As the hot fluid expands adiabatically, its thermal energy is converted 
to bulk kinetic energy. 

Since all the above authors assume an ultrarelativistic equation of state for the fluid, the 
final Lorentz factor their solutions predict for the fluid is formally infinite as the fluid never 
cools. They avoid this difficulty by following fluid elements in the flow only to the point 
where the fluid temperature becomes nonrelativistic and approximating the final coasting 
Lorentz factor as the one given by their solutions at that point. They thus find that the final 
Lorentz factor of a given fluid element scales as 7q +v ^ where 70 is the Lorentz factor acquired 
by the fluid when it is shocked. This method cannot accurately account for acceleration that 
occurs around the time when the fluid cools to nonrelativistic temperatures and can only 
produce approximate relations for the energy and velocity of the ejecta. While Kikuchi 
and Shigeyama (2007) relax the assumption of an ultrarelativistic fluid in their work on 
this problem, they cannot completely characterize the acceleration while the fluid is cooling 
either. 

We approach this problem by introducing a new kind of self-similar solution for the 
cooling and expanding fluid. In this solution, we require that the fluid move at relativistic 
speeds but relax the assumption that the fluid be hot. We place the characteristic position at 
the point where the fluid temperature is transrelativistic. We thus exploit the self-similarity 
of the transition between hot and cold fluid in the flow rather than the self-similarity in the 
acceleration of the hot fluid. Indeed, this flow when taken in its entirety is not self-similar: 
the size scales that characterize the acceleration and the hot/cold transition evolve with 
time according to different power laws. In other words, the entire flow is a composite of 
two distinct self-similar solutions. In §2 we summarize the solution for the hot fluid, which 
gives the initial conditions for this new solution. In §3 we derive the new solution, and in 
§4 we describe the behavior of fluid elements in the composite solution. In §5 we explain 
changes in the flow's behavior for very shallow initial density profiles in the stellar envelope. 
In §6 we discuss the behavior of the flow at late times and relate the elements' final Lorentz 
factors and the Lorentz factors to which they were initially shocked. In §7 we find regions 
of the flow where spherical corrections are important, and in §8 we summarize our findings 
and discuss them in the context of previous work. We take the speed of light to be c = 1 
throughout our discussion. 
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2. Initial conditions: behavior of the hot fluid 

We are interested in the behavior at late times of a fluid flow which begins as a relativistic 
shock propagating through the outer layers of a star with a polytropic envelope. As long 
as the distance between the front of the flow and the original location of the star's surface 
is small compared to the star's radius, the geometry is planar. So we seek a self-similar 
solution to the following hydrodynamic equations representing energy, momentum, and mass 
conservation: 

^ t b 2 (e + (3 2 p)}+^[l 2 (3(e + p)}=0 (1) 
|[7 2 /3(e+p)]+^[ 7 2 (/3 2 e + p)]=0 (2) 

| M + jL (7/3n) = . (3) 

The solution we seek must be connected to the hot flow, whose behavior is well understood: 
Sari (2006) and Nakayama and Shigeyama (2005) derive the self-similar solution before the 
shock breaks out of the star and Pan and Sari (2006) derives the post-breakout solution. 
Here we simply state these results. We take R(t) to be the characteristic position in the 
solution, and we set R = and t — at breakout. We take T, P, and N to be the 
characteristic Lorentz factor, pressure, and number density; R and T are related in that 
R = ^\ - l/r 2 -1-1/ (2r 2 ). We define 

tt m tP tN m , 

— = , — = — m — k , — — = k (4) 

r 2 ' P N 2 v ; 

Here k gives the unshocked density profile in the stellar envelope: assuming gravity is con- 
stant in the star's outermost layers and unshocked pressure and density are related by a 
power law, the density is given by a power law p oc \x\~ k where x is the position relative to 
the star's surface. This implies — 3 < k < — 3/2 for degenerate and convective envelopes and 
k = —17/13 for Kramers opacity envelopes. We consider here the regime k<-(l + V3) /3, 
which includes all of these profiles; we explain this choice of maximum k in §3. We write the 
solutions in the form 

1 \x ) t) = -T*(t)g{ X ) (5) 
p(x,t) = P(t)f(x) (6) 



where the similarity variable is 



X = l + 2(m + l)-^ . (8) 
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Note that this expression is equivalent to 



X 



t — x 
t-R 



(9) 



taken in the limit where r > 1, or where t ~ R(l + 1/ (2(m + l)r 2 ), so that the flow's 
characteristic length scale is t — R = R/ (2(m + 1)T 2 ). 

Mathematically, the pre- and post-breakout solutions differ only in the ranges in x which 
apply. They are 



-oo < x < 1 
oo > x > , 



, -oo < gx < 1 

oo > gx > (gx) 



t < (pre-breakout) 
t > (post-breakout) 



(10) 
(11) 



where 

(g X ) = 4 + 2V3-2kV3 • (12) 

These ranges differ because the sign of R and the interpretation of the characteristic values 
T, P, N, and R change at breakout. Before breakout, R < and T, P, N, and R are 
associated with the shock front; after breakout, R > and T, P, N, and R are associated 
with a fluid element which has expanded by a factor of order unity. 

The solutions are then completely specified by an expression for m in terms of k, which 
gives the time evolution of the flow, and expressions for g, f, and h in terms of Xi which give 
the spatial profiles for the hydrodynamic variables: 



m 



(3 - 2^) 



k 



9 = 
f = 

-gx + (gx), 



gx - (9x\ 



-i + (gx) 
gx - (gx) 



(gx) 



(3-2\/3)fc 



-2V3)fc 



(2\/3-3)(2fc-l)fc 



\9X 



2 -l+ky/3-V3 



(13) 



(14) 



(15) 



(16) 



i - (gx) 

Thus defined, Eqs. 5, 6, 7 satisfy the hydrodynamic equations with the equation of state 
p = e/3. They are an accurate description of the flow with k < — (l + \^3) /3 only where 
the fluid is hot, or where p/n ^> 1. As the fluid expands and accelerates after breakout, 
it cools adiabatically from the back of the flow towards the front (Pan and Sari 2006). As 
a result, the above post-breakout solution holds only for a region at the front of the flow, 
and this region shrinks with time. This solution sets the boundary conditions for the new 
solution we seek: as we approach the vacuum interface at x — 0, where the fluid is still hot, 
the new and the old solutions must coincide. 
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3. Self-similar solution for the cooling fluid 



We are interested in the behavior at late times of a fluid flow that begins as a relativistic 
shock propagating through the outer layers of a star with a polytropic envelope. We under- 
stand the behavior of the part of the flow that is hot {p/n 3> 1) and therefore obeys the 
equation of state p = e/3: it follows the self-similar solution given in Eqs. 4, 10-16. We note 
for convenience in our discussion below that the self-similar variable \i as given by Eq. 9, is 
equivalent to 



taken in the limit where T ^> 1, or where t ~ R(l + 1/ (2(m+l)r 2 ). The implied characteristic 
length scale is t - R = R/(2(m + l)r 2 ). 

As the fluid expands and accelerates after breakout, it cools adiabatically from the back 
of the flow towards the front. At late times, then, the post-breakout solution of §2, which 
we will refer to here as the "hot solution," holds only for a region at the very front of the 
flow, and this region shrinks with time. The hot solution sets the boundary conditions for 
the new solution we seek: as we approach the vacuum interface at x — 0, the two solutions 
must coincide. 

In the new solution, which we will refer to as the "cooling solution," we must include 
cold fluid. We therefore use the equation of state 



rather than the ultrarelativistic p = e/3. Although this equation of state applies only to 
fluids with adiabatic index 4/3, our analysis can be easily modified to accommodate an 
arbitrary equation of state of the form 



where F is an invertible function since the hydrodynamic equations can be written in self- 
similar form with any such equation of state. 

We must also specify a characteristic scale and define the characteristic Lorentz factor, 
pressure, and number density to be consistent with this scale. We seek the profiles of the 
hydrodynamic variables in the region where the fluid temperature transitions from hot to 
cold; the natural scale for this transition is the distance 5 between the vacuum interface, 
where the fluid is hottest, and the point where the fluid temperature becomes nonrelativistic. 
We set this point to be where p/n = 1. Then the similarity variable is 



V = g(e-™) 



(18) 




(19) 



t — X 



(20) 



5 
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and, by analogy with the hot solution, we express 7, p, and n as 

T 2 (x,t) = ±T a (t)m (21) 

p(x,t) = P(t)f(0 (22) 

n(x,t) = N(t)0^ . (23) 

We choose T, P, and N, the new characteristic values of the Lorentz factor, pressure, and 
number density, to match the 7, p, and n values given by the hot solution where p/n = 1. 
We take Xcoid to be the value of the old similarity variable x corresponding to p/n = 1 in 
the hot solution: 

r 2 = r 2 s( Xcold ) (24) 

P = ^/(Xcoid) (25) 
N = N f/ 2 f °' d) . . (26) 
This choice of characteristic values dictates 

P = N . (27) 
In the limit of late times, when 5 <C R/T 2 and gx — (gx)o ^ (fl'X)o' -^ ( l s - ^4 and 17 give 

x_ (^)o 1 _ fo , 1 



2 + V/3W . (29) 



2(m + i)r 2 v J r 2 

Note that the characteristic scale t — R~ R/T 2 in the post-breakout solution for the hot fluid 
is different from the new scale 5. R is the location of a fluid element that has expanded by a 
factor of order unity since breakout. Because R evolves according to the finite characteristic 
Lorentz factor T, R lags farther and farther behind the front of the flow, where the Lorentz 
factors are arbitrarily large. Since T decreases with time as per Eq. 4, t — R increases with 
time. In the limit of late times, then, R lags far behind the portion of the flow where the 
fluid remains hot and t — R becomes much larger than the space occupied by the hot fluid. 
In other words, the scale t — R that characterizes the hot solution becomes irrelevant to the 
transition between hot and cold fluid that is of interest here. 

To get r and 5 as functions of time, we apply Eqs. 14-16 at the point p/n = 1. We 
use P/N = r/(3\/2), a relation that follows from the shock jump conditions applied in the 
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pre-breakout solution. 



1 = p = r Vff(xcoid)/(xcoid) , 3Q x 

n 3^/2 ^(Xcoid) 

r 1 f 3m i ,„ i m(2fc-l) \ -fc 

= ^= [<?(Xcoid)H 2 [ 5 ( Xcold ) . Xcold _ 2J-1-V3-+W3- (31) 



r 2 5(xcoid) 



= 



3^ 



(#(Xcold) 'Xcold -2)-WsW3 



m /( m fc ) 



(32) 
(33) 



This gives 



tr 
T 



J 



l + V3 + 3k 
1 + + (3 - 2VS) k 



(34) 



(35) 



1 + V3 + 3k 

at late times, when #(Xcoid) • Xcoid ^ (gx)o- Solving Eq. 30 for /(Xcoid) or /i(Xcoid)/ y/g(x<x>id) 
similarly gives 



b- — - — 



4k 



r- ~ ■ (36) 

1 + V/3 + 3A; 

Note that k — — (l + VSj /3 ~ —0.91 makes a and 6 diverge in Eqs. 34 and 36. This k 
marks a qualitative change in the behavior of the flow which we discuss further in § 5. 

We now proceed to solve the hydrodynamic equations for the relevant range in k. We 
use the equation of state Eq. 18 to rewrite Eq. 1, Eq. 3, and the difference equation obtained 
by subtracting Eq. 2 from Eq. 1. We take the limit 7> 1. We rewrite the differentiation 
operators as 



d =, d ■ d • d 1 

— = rw + p — h TV h - 

dt dT OP ON 5 
d 1 d 
dx 5 <9£ 

and substitute these and Eqs. 24-27, 34, and 36 to get 



(37) 
(38) 



= 6 2/ 



h 

~T/2 



t 



5T 2 



-g 



J 3 h 

4^r + -- 



g 2 



4/' 

+ — + 
g g 



h' 



3/2 



? s 



g' h nS h' 
+ 2/ + 



2 g s / 2 



1/2 



(39) 
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= (2a + b)[2gf + 



hg 



1/2 



+ 



5T 2 



h 



4 g 3 / 2 



+ f' + 



h! 



2gV 2 



? s 



2/ + 



h 



= ( a + b)h + 



Ag 1 / 2 



t 



5T 2 



,h h' 
-9^ + ~ 



(40) 



(41) 



r g. 

We substitute Eqs. 29 and 35 into Eqs. 39, 40, and 41 and integrate this ODE system 
numerically to produce the solution shown in Figures 1,2, and 3. These figures also compare 
our solution to a one-dimensional numerical simulation of a relativistic planar shock wave 
which accelerates through and breaks out of a k = —3 medium and then expands and cools 
through the transrelativistic regime. The initial conditions used in the simulation represent 
the shock as an unresolved contact discontinuity with r = 5 at x = — 10~ 5 at starting time 
t — — 1.0 x 10~ 5 . The simulation grid contains 400 fluid elements spread logarithmically over 
two orders of magnitude in x; our calculation follows these elements until time t = 2.6 x 10 8 . 

We can check that the behavior of this solution at large £ — where the fluid is very cold 
and where the hot solution and cooling solution differ most — is physical. Consider a fluid 
element many distance scales 8 behind the vacuum interface at position t — x ^> 5. This 
fluid element must have become cold at some time t co id «C t; as a result, it has long since 
stopped accelerating and has spent most of the time interval t — t co id coasting at its current 
Lorentz factor 7. Then this fluid element has 

t-(t- tcoid) V 1 - V7 2 */ (f 2 9) 2 - Vs 



t 



x 



t/t 2 ■ (2 + V3) g 



(42) 



V3 



(43) 



so at large £ we expect 

m = 2 

We cannot get exact relations for / and h in the large £ limit in this way because p and n 
change significantly while the fluid element finishes its acceleration. However, we can check 
the scalings of / and h with £. Because the fluid elements far behind the front are coasting 
with Lorentz factors that are virtually constant in time, the volume of each fluid element 
increases linearly with time. This implies 



n oc t 



p oc n 4 / 3 oc t 



f-4/3 



(44) 



for a single fluid element. From the definitions of T, P, and N in Eqs. 24, 25, and 26, we 
know 7(£ = l)/r, p(£ = 1)/P, and n(£ — 1)/N are constant in time. Then for a single fluid 
element, 



V ~ Mucoid ) 



^-cold 



-4/3 



oc i 



6+4/3 

cold 



(45) 
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log 10 (p) 




iog 10 © 



Fig. 1. — Profile of the pressure p as a function of the similarity variable £ for k = —3. The 
dashed line is the hot solution valid for the fluid near the front, at small £; the solid line 
is the cooling solution. Data from numerical simulations are shown as crosses. In order to 
cover a substantial range in £, data from six p vs. £ profiles corresponding to different times 
in the same simulation run are shown. The data agree well with the cooling solution. The 
"tails" at the ends of the numerical simulation profile data are due to edge effects at the ends 
of the simulation grid that are not self-similar. The overall y-axis normalization is arbitrary, 
but the relative normalizations of the hot solution, the cooling solution, and the numerical 
simulations are correct. 
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log 10 (n) 




log 10 (f) 



Fig. 2. — Same as Figure 1 for the number density n rather than p. 
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log 10 (y) 




iog 10 (f) 



Fig. 3. — Same as Figure 1 for the Lorentz factor 7 rather than p. 
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Since 



9 = 



2 7 2 r 2 (t cold ) 



oc t 



2a 

cold ) 



(46) 



f 2 (t) f 2 (t) 



we have 



6+4/3 

p OC q 2a 



(47) 



A similar calculation yields 




(48) 



That the relations in Eqs. 43, 47, and 48 hold at large £ is shown in Figure 4. 



4. The composite solution 



As discussed in § 2 and 3, the hot fluid close to the vacuum interface and the cooling 
fluid further back in the flow obey two different self-similar solutions in which both the time 
evolution and the physical interpretation of the characteristic length scale, Lorentz factor, 
pressure, and density differ. In other words, we describe the entire post-breakout flow — 
which, when taken as a whole, is not self-similar — by a composite of two self-similar solutions 
built up around the two different length scales which characterize different portions of the 
flow. This to our knowledge is the first such 'composite solution' found. Chevalier (1982), for 
example, also uses two self-similar solutions in a description of a single flow — specifically, the 
interaction between an expanding shell of fluid and a stationary external medium. However, 
his solutions have the same characteristic length scale and time evolution; they differ only 
in the shapes of their profiles and in the disjoint regions of the flow in which they operate. 
They may be considered as a single self-similar solution in which the pressure, velocity, and 
density profiles are piecewise functions of position. 

We can check that the behaviors of fluid elements and sound waves in the hot solution 
and the cooling solution are consistent. The characteristic position R in the hot solution 
moves backwards relative to the vacuum interface with time because T decreases with time. 
By contrast, the characteristic position t — 5 in the cooling solution moves forwards relative 
to ct because fluid elements at the back of the solution cool faster than those at the front, 
and t — S marks the location of a fluid element that has just cooled. We confirm the forward 
motion by looking at Eq. 35, which indeed gives t5/5 < for the range of k of interest 
(k < — (l + y/3) /3). We expect fluid elements in the cooling solution to move backwards 
in the solution, or towards larger £: every fluid element must eventually finish accelerating 
and become cold, so the point t — 8 that marks the hot/cold transition must overtake every 
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Fig. 4. — Verification of the scalings of g, /, and h with £ at large £, or cold fluid temper- 
atures. The functions plotted (dashed line for g, dotted line for /, solid line for h) were 
obtained via numerical integration of the ODEs in Eqs. 39, 40, 41. They show that Eqs. 43, 
47, 48 — relations derived for fluid which has finished accelerating — are valid at large £. In 
particular, g£ approaches the expected value 2 — = 1CT - 572 . In this calculation we used 
k = -3. 
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fluid element. Indeed, the time derivative of £ following a fluid element, 

°L = (± + t fcII±\ ( =i(l^-*i\ (49) 
Dt \dt V 7 2 <W t\ U 6 J • { ' 

is always positive since t5/5 < 0. Fluid elements in the hot solution move forwards with 
time, towards smaller gx, since they accelerate while R decelerates. The advective time 
derivative of Xi 

is always negative since gx > (gx)o > 4 + 2\^3 everywhere. 



If the proper sound speed in the fluid is (3 S = ^ A/3f 1 l 2 g 1 l i h 1 ^ 2 , then the motion of a 
sound wave in the cooling solution is given by 



dln£± _ t 1 
dint £±5 



( dx ± A 2- y/3 f lTPs \ t8 



where the signs denote forward- and backward-propagating sound waves. Again, tS/S < 0, 
so d£±/dt > everywhere and all sound waves move backwards in the cooling solution. In 
other words, all fluid elements are disconnected from the vacuum interface. In the sense that 
the front is disconnected from the fluid far back in the flow, the cooling solution is similar 
to Type II solutions. However, in contrast to the usual Type II scenario, there is no sonic 
point constraining the solution. 

There is a caveat in this composite view of the flow for the density profiles k < 
— (l + v^3) /3 discussed in the preceding sections. At early times when T ^> 1, all the 
fluid elements which have accelerated by at least a factor of order unity are those with 
p/n>T; slower and cooler fluid elements still have very nearly the same temperatures and 
speeds as they had just after being shocked. So the only fluid elements with p/n of order 
unity or smaller are those which were never shocked to relativistic temperatures, and the 
part of the flow moving relativistically can be described with the hot solution alone. At 
late times when the cold solution is relevant, the hot solution applies only in a region at the 
front of the flow whose size is much smaller than the characteristic scale R/T 2 . In this small 
region, the profiles of the hydrodynamic variables in hot solution appear very nearly scale 
free, and the cooling solution approaches the resulting power-law profiles toward the front 
of the flow, in the limit of small £. Then the hot solution is unnecessary to the description 
of the flow at these late times, so we can think of the flow as following the hot solution at 
early times and transitioning to the cooling solution when T = 1. Figure 5 is a schematic of 
the composite flow at these late times. 
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lo 9io (gx - (gx)o) 



o 



Fig. 5. — Pressure profiles for initial density profile k = —3 at late times when the cooling 
solution applies. The dashed curve is the hot solution and the solid curve is the cooling 
solution. The front of the flow, where gx — (gx)oi is toward the left. As in Figure 1, 
the overall |/-axis normalization is arbitrary but the relative normalizations of the hot and 
cooling solutions are correct. The 'kinks' in the hot and cooling solution curves occur at 
their respective characteristic scales t — R and 5. Fluid far back in the flow, where £ 1 or, 
equivalently, gx — (<?x) ^ (9X) co \d ~~ (flX)o> follows the cooling solution. This cold fluid falls 
behind the position predicted for it in the hot solution because it is no longer accelerating. 
As a result, the hot solution gives artificially high pressures for fluid at the back of the 
flow. The hot solution applies only at the front of the flow in a region small compared to 
both characteristic scales. As this region corresponds to the limit in which we set the cooling 
solution to match the hot solution, the entire flow may be described with the cooling solution 
alone at these times when T < 1. 
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However, a composite of the hot and cooling solutions is essential in describing the flow 
for density profiles — (l + /3 < k < 0, which we discuss in the next section. 

5. Solutions when - (l + v 7 ^) /3 < k < 

We restricted our discussion of the hot and cooling solutions in § 2, 3, 4 to initial 
density profiles k < — (l + \/S) /3. That the initial density decrease toward the star's 
surface requires k < 0, but so far we have neglected the interval — (l + a/3) /3 < k < 0. To 
describe the flow for these k, we return to the singularity at k — — (l + V3) /3 noted in § 3 
and examine the behavior of (g-x) cold = #(Xcoid) • Xcoid- 

Eq. 30 implies 

dln (ffxLid = (^x)coid ~ (9X)o (ffX)coid ~ 2 / 52 n 

d^t (g X ) cold ( 5X ) cold -2(3 + 4/V3) 

For — (l + \/3) /3 < < 0, integration of Eq. 52 shows that (gx) co \d * s a double-valued 
function of time until the entire flow cools. As t — > from above, (gx) C oid — * {(fi , x)o'°°}; 
for positive times, the smaller value of {gx) co id increases and the larger value decreases with 
time until (<7x) C oid becomes single- valued at (<?x) co id = 2 (3 + 4/ \/3) , which occurs at a finite 
time. In other words, fluid cools quickly both at the front, where the expansion timescale is 
shortest, and at the very back, where the shocked fluid is coldest. 

Since (gx) C oid becomes single- valued at a (<7x)coid of order unity, the value of g corre- 
sponding to the larger value of (<7x) co id * s al wavs of order unity or smaller. But g — 1 for a 
fluid element which has accelerated by a factor of order unity, so this larger (gx) co id tracks 
fluid elements which cool before or just as they manage to accelerate by a factor of order 
unity. These fluid elements cannot have been shocked to relativistic temperatures, so we 
focus on the evolution of the smaller value of (gx) co \d as ^ moves through the front of the 
solution. In contrast to the k < — (l + VS) /3 case, the fluid near the front cools from the 
vacuum interface backwards, from the fastest-moving toward slower-moving fluid, and the 
fluid farthest forward is coldest: T decreases with time and P, N, 5 increase with time. 
Fluid elements just behind the front — specifically, those with gx — (gx)o ^ 2 + 2 /a/3 + 2k\/3 
— always have (<7x) cold — (<7x)o> so they obey Eqs. 34-36. Then the discussion of the cooling 
solution in § 3 applies in the limit of small positive times, when the front of the flow is cool- 
ing, except in that the initial conditions used to integrate Eqs. 39-41 are set by matching to 
the hot solution far back in the cooling flow, at £ ^> 1, rather than at the vacuum interface. 
Between breakout and the time when all the fluid cools, the flow contains both hot and cold 
fluid moving at relativistic bulk speeds; it therefore follows a true composite of the hot and 
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cooling solutions. Figure 6 shows a schematic of this composite solution. 

As in the k < — (l + \/3) /3 case, fluid elements move forwards in the hot solution 
with time according to Eq. 50. No qualitative change in the time evolution of the scale R 
—or, therefore, in the behavior of fluid in the hot solution — occurs as k increases through 
— (l + \/3) /3. In contrast to the k < — (l + Vs) /3 case, fluid elements move forwards 
in the cooling solution as well. Eq. 42, which still holds in the limit t — x <C 5, implies 
K > 2 - y/3, while Eq. 35 implies t8/8>\. Then from Eq. 49 we have D£/Dt < 0. 

For k — — {l + y/3) /3, (gx) = 3 (2 + 4/>/3), so (gx) co id is single- valued and (gx) coM - 
2 oc i -1 . Then there is a finite time at which (<7x) cold reaches (<?x) , 5 shrinks to 0, and 
the entire flow cools. Also, at any t > 0, Eq. 15, 16 imply that the temperature at the 
front of the flow, where gx — (gx)o ^ (fl/Oo' is roughly independent of position, and the 
temperature at the vacuum interface is finite; all the fluid at the front cools at the same 
rate. Though Eq. 35 does not work for this k, we know S decreases monotonically to with 
time, so Eq. 49 gives D^/Dt > 0. So fluid elements move backwards in the cooling solution 
until (^x)coid = (9X) and 5 = 0. 



6. Behavior of fluid elements at late times 

Earlier analytic work has established that the final Lorentz factor 7fi na i of a given fluid 
element should scale according to 

7final = Kt££ d (53) 

where 7 s hocked is the fluid element's Lorentz factor immediately after it is shocked in the 
pre-breakout flow and the coefficient K is independent of 7 s hockcd (Johnson and McKee 1971; 
Pan and Sari 2006). Tan et al. (2001) have found numerically that K ~ 2.6 for k = —3. 
They note, and we confirm from our own experience, that it is difficult to continue numerical 
simulations until the very end of the fluid acceleration since the conversion of thermal to 
bulk kinetic energy is quite slow: Figure 4 indicates significant acceleration until £ ~ 10 10 . 
To estimate the coefficient, Tan et al. (2001) applied correction factors to their simulation 
results of up to ~50% for fluid elements with final Lorentz factors of order ~ 10 3 . 

We can find 7fi na i for a given fluid element directly from our pre- and post-breakout 
solutions. To track the acceleration of the fluid element while it is hot, we take the advective 
time derivative of 7 in the pre- and post-breakout solutions for the hot fluid and integrate 
with the proper limits. 

^7_7 (VS~S)k 
Dt t gx -A-2VS 
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Fig. 6. — Same as Figure 5, but for initial density profile k = —3/4 at early times, before 
(<?x)coid becomes single- valued. Here the cooling solution applies at the front of the flow; the 
hot solution applies further back in the flow. Because the cold fluid has stopped accelerating, 
the hot solution gives artificially high pressures for fluid towards the front of the flow where 
£ <C 1 or gx — (<?x)o ^ (fl'X)coid — (fl'X)o- However, the cooling solution misses the transition 
between regions of hot fluid which have and have not accelerated significantly since being 
shocked; this transition occurs at x ~ 1, or gx—(9X)o °f or der unity. An accurate description 
of the entire flow requires both the hot and cooling solutions. 
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Dgx = 1 (2 - gx) (_9X ~ (gx)g) 
L>7 (v^-3) fc 

£>#x 7 (2 -gx)(gx- (gx) ) 

Before breakout, the fluid element's gx goes from gx = I when it 
at breakout. After breakout, the fluid element's gx goes from gx - 
when the fluid becomes cold, we have 

(V»-3)fc / ny — 2 \ (9X) -2 

7 = 7 , ocM ( ta)o -l) M -=(^-) . (57, 

To find the correct at which to evaluate the above, we set C = ^f/h to be the temperature 
(up to a factor of 3) at the gx of interest and use Eqs. 15, 16 to express {gx — 2)/(<7X — (gx)o) 
in terms of C . This gives 

7 = C-^7i2£d • (58) 

To this we add the extra factor given by the cooling solution to get the coefficients K shown 
in Figure 7. In particular, 

Tfinai = l-967i£Ld k = - 3 (59) 

7^1 = 2.717^ k = -3/2 . (60) 

This result is close to the results of Tan et al. (2001), who find a coefficient of ~ 2.6 when 
k = —3. Note that p/n = 1 corresponds to C _v ^ = 0.149: 7 grows by a factor of ~ 15 after 
the fluid element becomes nominally cold. 

The growth of 7 as a function of the temperature for a single fluid element is shown in 
Figure 8, which also shows good agreement between the cooling solution and direct numerical 
simulations of the hydrodynamic equations. Because the Lorentz factors near the front of 
the flow in particular become very large at late times, it is difficult to produce numerical 
simulations that remain accurate as the fluid cools all the way to p/n <^ 1. As a result, the 
numerical simulation shown in Figure 8 cuts off while the fluid Lorentz factor is 9% smaller 
than the final Lorentz factor predicted by the cooling solution. 



(55) 
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is shocked to gx — > —00 
-> 00 to gx ^ (gx) co i d - So 



7. Effects of spherical geometry 

We can estimate the ranges of initial positions and Lorentz factors for which corrections 
to our planar solutions due to the star's spherical geometry are important. Spherical geom- 
etry significantly affects a given fluid element's acceleration if the distance between the fluid 




Fig. 7. — Coefficient K for the final Lorentz factor as defined in Eq. 53 plotted as a function 
of density profile k. The crosses were computed from the cooling solution; the line connecting 
them is included to guide the eye. 
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Fig. 8. — Lorentz factor 7 of a single fluid element as a function of the temperature p/n 
of that fluid element. Lower temperatures and later times are towards the left. The solid 
line is the self-similar solution; the points are the results of numerical simulations. Both 
calculations were done for k = — 3. As the fluid element becomes cold, the evolution of its 7 
with p/n deviates from the power law seen at high temperatures. 
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element and the star's center doubles before the fluid element finishes accelerating. Once the 
fluid element has traveled a distance comparable to i?*, the star's radius, it has expanded 
significantly in directions perpendicular to its motion; planar solutions cannot account for 
this transverse expansion. For our estimate we therefore check which fluid elements have 
cooled to p/n ~ 1 by time t ~ R*. 

Consider a fluid element whose position before being shocked is rr init < 0. According to 
the relativistic shock jump conditions, immediately after this fluid element is shocked it has 

Tinit ~ T ~ ( X -^\ m ' 2 (61) 
mit \ X re \ I 



(") 

\n/ ini 



where Y is the shock Lorentz factor in the pre-breakout solution and x m \ is the position of 
the shock when it first becomes relativistic, or when T ~ \/2. Right after breakout our fluid 
element follows the hot solution; when it reaches x = 1, it has accelerated and expanded by 
a factor of order unity since being shocked, so it still has p/n ~ (^init/^rei) - "^ 2 - Also, at this 
point T ~ 7i n it so the time is~ |ximt|- We track the temperature of the fluid element over 
time until it reaches p/n = 1 by taking the advective time derivative of p/n: 

DHpN = k(V3-l) ^ y/3-3 

Dint 3 X _ 4 _ 2v /3 6 1 ' 

where the final equality holds while the fluid element is at gx — (gx)o ^ (fl'X)o — that is, 
near the front of the flow, where most of the acceleration takes place. We then have 

(* H E) ( « . m 

\nJ Vn/init v|x in it|/ 
We now impose the condition p/n ~ 1 at or before t ~ i?* to get the scaling 

, . 1+ 1+V3 

^< • (64) 

Sphericity corrections are unimportant for fluid elements which satisfy this condition. For 
> k > — (l + VSj /3, the exponent on the right hand side is not positive, so spherical 
corrections are unimportant as long as |x rc i| < i?*. On the other hand, for k < — (l + y/3) /3 
spherical corrections are important for a layer of fluid initially adjacent to the star's surface: 
in our self-similar solutions for these k values, fluid elements that start arbitrarily close to 
the star's surface will take arbitrarily long to cool. 

The equivalent condition on 7i n i t , the Lorentz factor of a given fluid element just after 
being shocked, is 

*£rel 



l+%/3+3fc 

< (Tinit) (3+V5)fe ; (65) 
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the one for 7fi na i, the fluid element's Lorentz factor after it finishes accelerating, is 

|x r el| , i+Vs±3k 

' < (Tfinal) ^ • (66) 



R 



8. Summary and discussion 

We have derived a new self-similar solution, the cooling solution, for the flow that results 
when a relativistic shock breaks out of a polytropic envelope. The cooling solution is based 
on our identification of the characteristic position with the point where the fluid cools to 
nonrelativistic temperatures. The cooling solution shows that the transition between hot 
and cold fluid in the flow is self-similar even though this transition is not included — indeed, 
is not self-similar — in the old post-breakout solution for the hot fluid alone. We present a 
description of a non-self-similar flow using a composite of two distinct self-similar solutions 
in which the time evolution and physical interpretation of the characteristic scales differ. As 
the envelope's original density profile k increases past k — — (l + \/3) /3, the flow dictated 
by the composite solution changes qualitatively from one where the slowest fluid at the back 
of the flow cools most quickly to one where the fastest fluid at the front of the flow cools 
most quickly. 

The cooling solution allows accurate calculation of the final Lorentz factors of the 
shocked fluid elements. Given a stellar model for a core-collapse supernova progenitor and 
an input explosion energy, we can extract the initial density profile of the progenitor enve- 
lope and the initial shock velocity and use the cooling solution to find the Lorentz factor 
and kinetic energy profiles in the relativistic ejecta after the ejecta finish accelerating. This 
provides an accurate value for the energy available to produce observable lightcurves via 
interaction between the ejecta from the model explosion and the progenitor's surroundings. 

Kikuchi and Shigeyama (2007) also investigate cooling in the flow produced after a 
relativistic shock breakout. They focus on a k = —3 density profile and use as starting point 
the work of Nakayama and Shigeyama (2005), who found a numerical self-similar solution for 
the hot planar flow in a Lagrangian framework. While Kikuchi and Shigeyama (2007) also 
describe the cooling flow with a system of ODEs, they conclude that no self-similar solution 
exists for the cooling flow: instead of adopting a new characteristic scale in writing the ODEs 
for the cooling flow, they retain the characteristic scale relevant to the hot solution. Because 
they do not recognize the self-similarity of the flow, and because they do not integrate until 
the fluid has cooled enough, they find that initially hot fluid does not stop accelerating. As 
a result they cannot express the final Lorentz factors of given fluid elements in terms of the 
initial ones. 
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In addition, Kikuchi and Shigeyama (2007) focus on corrections to their planar flow due 
to spherical geometry. They perturb the hydrodynamical equations to lowest order in a new 
variable equal to the position coordinate scaled to the stellar radius. Their non-self-similar 
solutions to the perturbed equations agree well with their numerical simulations only at very 
early times, before the fluid has expanded by a factor of ~2. They look at sphericity effects 
in the cooling flow using numerical simulations and find that for fluid whose initial fractional 
distances from the star's surface to its center lie between about 4 x 10~ 4 and 3 x 10~ 3 , 
sphericity effects change the final velocities by factors of order unity in simulations with 
k = —3, T = 10 5 and 10 6 at t — 1, and i?* = 3.3 in units where c = 1. Although such a fast 
shock is unrealistic since it would have T ~ 10 4 or 10 5 at the star's center, our scaling in §7 
also indicates that sphericity corrections would be important for these parameters. 

Tan et al. (2001) also consider sphericity corrections for the density profile k = —3, but 
they treat the opposite limit of mildly relativistic shocks with T — 1 ranging from about 
7 x 10~ 3 to 0.28 when the shock is halfway from the star's center to its surface. Only 
their two most energetic numerical calculations of sphericity effects attain shock Lorentz 
factors T > \/2~, so nonrelativistic estimates of sphericity corrections are relevant for most of 
the regime they consider. However, those two most energetic shock calculations show that 
for the relativistic fluid elements shocked to initial Lorentz factors between about y/2 and 
4, sphericity decreases the final velocities by twenty to thirty percent. Our scaling in §7 
indicates that for those two scenarios, sphericity corrections should be important for fluid 
elements shocked to Lorentz factors larger than about 1.7 and 2.4 respectively; this is roughly 
consistent with the findings of Tan et al. (2001). 
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